% Natasha Savage

lp = 0;
react_t = 0;
et = 0.0001;

    while  lp ~= 2
        %Masha dev sept 2012 Added RecCmplxic,RecCmplxn,RecAn species to all Euler function calls 
        [Cdc42Th, BemGEF42h, BemGEFh, BemGEFch, Cdc42Dh, GDI42h, GDI42ch, GDIch, Cdc42Dich, GDI42ich, RecCmplxich,RecCmplxh,RecAh] = ...
            Euler_step_GDI42ic(Cdc42T,  BemGEF42,  BemGEF,  BemGEFc,  Cdc42D,  GDI42,  GDI42c,  GDIc, Cdc42Dic, GDI42ic,RecCmplxic,RecCmplx,RecA,    dt2);
        
        [Cdc42Th2,BemGEF42h2,BemGEFh2,BemGEFch2,Cdc42Dh2,GDI42h2,GDI42ch2,GDIch2, Cdc42Dich2,GDI42ich2,RecCmplxich2,RecCmplxh2,RecAh2] = ...
            Euler_step_GDI42ic(Cdc42T,  BemGEF42,  BemGEF,  BemGEFc,  Cdc42D,  GDI42,  GDI42c,  GDIc, Cdc42Dic, GDI42ic,RecCmplxic,RecCmplx,RecA,    dt2/2);
        [Cdc42Th2,BemGEF42h2,BemGEFh2,BemGEFch2,Cdc42Dh2,GDI42h2,GDI42ch2,GDIch2, Cdc42Dich2,GDI42ich2,RecCmplxich2,RecCmplxh2,RecAh2] = ...
            Euler_step_GDI42ic(Cdc42Th2,BemGEF42h2,BemGEFh2,BemGEFch2,Cdc42Dh2,GDI42h2,GDI42ch2,GDIch2,Cdc42Dich2,GDI42ich2,RecCmplxich2,RecCmplxh2,RecAh2,dt2/2);
        
        Cdc42Te    = max(max(abs(Cdc42Th   - Cdc42Th2   )));
        BemGEF42e  = max(max(abs(BemGEF42h - BemGEF42h2 )));
        BemGEFe    = max(max(abs(BemGEFh   - BemGEFh2   )));
        BemGEFce   = max(max(abs(BemGEFch  - BemGEFch2  )));
        Cdc42De    = max(max(abs(Cdc42Dh   - Cdc42Dh2   )));
        GDI42e     = max(max(abs(GDI42h    - GDI42h2    )));
        GDI42ce    = max(max(abs(GDI42ch   - GDI42ch2   )));
        GDIce      = max(max(abs(GDIch     - GDIch2     )));
        Cdc42Dice  =        abs(Cdc42Dich  - Cdc42Dich2 )  ;
        GDI42ice   =        abs(GDI42ich   - GDI42ich2  )  ;
        RecCmplxice=        abs(RecCmplxich-RecCmplxich2)  ; %Masha dev sept2012
        RecCmplxe  =max(max(abs(RecCmplxh  -  RecCmplxh2))); %Masha dev sept2012
        RecAe      =max(max(abs(RecAh      -      RecAh2))); %Masha dev sept2012
        
        r = max([Cdc42Te,BemGEF42e,BemGEFe,BemGEFce,Cdc42De,GDI42e,GDI42ce,GDIce,Cdc42Dice,GDI42ice,RecCmplxice,RecCmplxe,RecAe]);

        % doing the step again with different stepsize        
        if r > et && lp ~= 2
            dt2 = dt2/1.5;     

            if lp == 1
                lp = 0;
            end
            % moving on
        else

            Cdc42T   = Cdc42Th2;
            BemGEF42 = BemGEF42h2;
            BemGEF   = BemGEFh2;
            BemGEFc  = BemGEFch2;
            Cdc42D   = Cdc42Dh2;
            GDI42    = GDI42h2;
            GDI42c   = GDI42ch2;
            GDIc     = GDIch2;
            Cdc42Dic = Cdc42Dich2;
            GDI42ic  = GDI42ich2;
            RecCmplxic= RecCmplxich2;
            RecCmplx = RecCmplxh2;
            RecA     = RecAh2;
            
            % adding the step on 
            react_t = react_t + dt2;
            % setting new stepsize for next step
            dt2 = 1.5*dt2;
            
            if lp == 1
                lp = 2;
            end
            
        end
            
        % Making sure we stay within the dt interval
        if react_t + dt2 > dt && lp ~= 2
            dt2 = dt - react_t;
            lp = 1;
        end            
            
        
        %checking for negative values and reducing tolerance 
        if lp == 2

           if (~isempty(find(Cdc42T<0))  | ~isempty(find(BemGEF42<0)) | ~isempty(find(BemGEF<0)) ...
             | ~isempty(find(BemGEFc<0)) | ~isempty(find(Cdc42D<0))   | ...
               ~isempty(find(GDI42<0))   | ~isempty(find(GDI42c<0))   | ~isempty(find(GDIc<0)) ...
             | ~isempty(find(GDI42ic<0)) | ~isempty(find(Cdc42Dic<0)) ...
             )
           
            lp = 0;
            react_t = 0;
            et = et/10;
           
            fprintf(neg_store,' %5.4f \t%5.4e \t%5.4e \n', curr_t,et*10,et);  
           end
           
        end    

    end
    


    %% diffusion in cargo step
    

    
    
  





